Defining the domain

First, we create a proper domain object exploiting the function create_mesh which takes a R named list as parameter. Such named list should contain geometrical information of the domain:

• nodes, a 2-columns matrix containing the coordinates of the nodes of the mesh.

• elements, a 3-columns matrix that for each row \(i\) contains the indexes of the mesh nodes that are vertices of the \(i\)-th element.

• boundary, a 1-column matrix that indicates the nodes on the boundary of the domain.

femR contains the unit_square named list that can be used to perform some tests.

data("unit_square", package="femR")
unit_square <- Mesh(unit_square)
plot(unit_square)

Moreover, femR package provides an utility that reads files storing such geometrical information and provided by third-party software such as FreeFem++. Indeed, read_mesh(filename) reads a .mesh file from FreeFem++ and returns a named list that can be passed to create_mesh function. The following chunk, that is not executed, shows how read_mesh works.

filename <- "path/to/domain.mesh" #domain.mesh stored by the FreeFem++ utility savemesh() 
domain <- read_mesh(filename)
mesh <- Mesh(domain)

Defining the solution of the PDE and the differential operator

First, we define initialize a FunctionSpace object passing the mesh tha we previously build and the finite element order. Note that, the package provide first order and second order finite elements. First order finite element is the default parameter. Then, we define a Function, solution of the PDEs, belonging to the function space Vh.

## solution of the PDE
Vh <- FunctionSpace(unit_square, fe_order=2)
f <- Function(Vh)

Then, we define the differential operator in its strong formulation.

L <- -laplace(f) # L <- -div(I*grad(f))
                 # In general:
                 # L <- (-mu)*laplace(f) + dot(b,grad(f)) + c*f

Note that, according to the well-known identity \(div(\nabla f) = \Delta f\), we could have written L<--div(I*grad(f)), where I is the 2-dimensional identity matrix. Indeed, femR package let you define every term of a generic second-order linear differential operator \(Lf= -div(K\nabla f) + \mathbf{b}\cdot \nabla f + a f\) as follows:

• Diffusion: -div(K*grad(f)) or (-mu)*laplacian(f), where either K or mu represents the diffusion matrix in an anisotropic problem or the diffusion coefficient in a isotropic problem, implements the first term of the differential operator according to the problem at hand.

• Advection: dot(b,grad(f)), where b is a two-dimensional vector, implements the second term of differential operator.

• Reaction: c*f, where c is a scalar, implements the last term of the differential operator.

Although, the interface let the user define a generic second-order linear differential operator: -div(K*grad(f)) + ..., at this stage the package has been tested only in the isotropic case, i.e K=I, where I is the two-dimensional identy matrix.

Building the PDE object

We define the forcing term of the differential problem and the Dirichlet boundary condition as a simple R functions. Such functions will be provided as parameters to the object which encapsulates the concept of PDEs.

## forcing term
u <- function(points){
    return(8.*pi^2* sin(2.*pi*points[,1])*sin(2.*pi*points[,2])) 
}

## Dirichlet BC
dirichletBC <- function(points){
  return(matrix(0,nrow=nrow(points), ncol=1))
}

Computing the FE solution

Finally, we build a pde object passing the differential operator L, as first parameter, the forcing term u, as second parameter, the Dirichlet boundary condition dirichletBC, as third parameter:

## create pde
pde <- Pde(L, u, dirichletBC)

We can compute the discrete solution of the Poisson problem calling the solve method:

## solve problem
pde$solve()

Since, we are dealing with a simple differential problem and we know the analytic form of the exact solution, we can compute the \(L^2\) norm of the error exploiting on the mass matrix provided by the pde object:

exact_solution <- function(points){
    return( sin(2.*pi*points[,1])*sin(2.*pi*points[,2]))
}
u_ex <- as.matrix(exact_solution(pde$get_dofs_coordinates()))
error.L2 <- sqrt(sum(pde$get_mass()%*%(u_ex-pde$solution())^2))
cat("L2 error = ", error.L2, "\n")
#> L2 error =  0.0004136347

Graphics

The R package plotly is exploited to plot the computed solution.

## plot solution 
plot(f)

# contour 
contour(f)

Note that, the base::plot function and graphics::contourfunction have both been overloaded and return a plotly object that can be modified as one wishes.

plot(f, colorscale="Jet") %>% layout(scene=list(aspectmode="cube"))